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SUMMARY 

Two numerical codes have been developed for the calculation of three-dimensional 
nozzle exhaust flow fields associated with hypersonic airbreathing aircraft. Both codes 
employ reference plane grid networks with respect to three coordinate systems. Pro- 
gram CHAR3D is a characteristic code utilizing a new wave preserving network within 
the reference planes, while program BIGMAC is a finite difference code utilizing conser- 
vation variables and a one-sided difference algorithm. Secondary waves are numerically 
captured by both codes, while the underexparision shock and plume boundary are treated 
discretely. The exhaust gas properties consist of hydrogen-air combustion product mix- 
tures in local chemical equilibrium. Nozzle contours are treated by a^ newly developed 
geometry package based on dual cubic splines. Results are presented for simple config- 
urations demonstrating two- and three-dimensional multiple wave interactions. 

INTRODUCTION 

Hypersonic^ air craft with airbreathing propulsion will require a high degree of 
engine/airframe integration in order to achieye optimized performance. The engine 
exhaust flow, because of physical area^limitations, will generally be under expanded at the 
nozzle exit, and in order to obtain maximum propulsive efficiency, the vehicle afterbody 
undersurface is used to provide additional expansion. This results in a three-dimensional 
nozzle flow whose boundaries are defined both by the solid boundary of the nozzle wall and . 
by the boundary separating the nozzle flow from the vehicle external flow. A typical 
exhaust nozzle (fig. l) may be characterized as having nozzle modules with cross sections 
which are rectangular in shape. These nozzles may be arranged in multiples and dis- 
charge into a common nozzle. The flow fields to be analyzed start at the combustor exit 
and each module may be analyzed individually until its merger with adjacent modules and 
the external flow field. 

*This research was performed under Contract No. NAS 1-12726. 


In developing a numerical model for this flow field, the following dominant features 
must be accounted for; 

(1) The flow properties at the combustor exit are highly nonuniform. Burning and 

mixing in the combustor yield regions of highly varying composition, temperature, and 
stagnation properties. In addition, shock waves are produced in the vicinity of the injec- .. 
tors. Although the strength of these waves decays rapidly as they propagate through the 
burner, they are generally present at the burner exit and must be accounted for. t 

(2) The exhaust gas mixture consists of hydrogen-air combustion products and sig- 
nificant burning may still occur in the initial regions of nozzle expansion. 

(3) The flow field geometry is quite complex. The engine modules consist of mul- 
tiple surfaces with sharp interior corners, and flow fences to contain the external exhaust 
flow may be present. 

(4) The interior nozzle flow field is dominated by complex wave interactions with 

waves generated and reflected off multiple surfaces. In addition, sharp interior corner ' ' 
regions must be accounted for. - 

(5) The nozzle exhaust flow interacts with the nonuniform vehicle external flpw 

field. This complex interaction for underexpanded exhaust flows results in an expansion . . 
system propagating toward the vehicle undersurface from the cowl trailing edge and a 
spanwise expansion generated by the sidewall interaction. An underexpansion shock prop- 
agates outward into the nonuniform vehicle external flow, and the exhaust and external 
flow are separated by a plume boundary. In addition, pressure and flow deflection mis- 
match between adjacent modules may occur, resulting in a spanwise multiple shock 
system. ' 

To best accommodate highly rotational variable composition flow fields, a grid nety - 
work which follows streamlines is preferred. For nonstreamline networks, large errors. . 
may be associated with streamline interpolation procedures for nonequilibrium flow calT ^ 
culations, as discussed by Sedney (ref; 1). For two-dimensional flow fields, a grid net- 
work following the flow streamlines is readily obtained. Such a system is employed in ; 
references 2 and 3 for the calculation of chemical reacting nozzle flow fiel^ and super- 
sonic combustor flow fields, employing a "viscous" characteristics technique. In this , v 
approach, a uniform marching step Ax is taken, new streamline grid pointy, are obtained,, 
and characteristic data are obtained by interpolation on the initial data line. Such a / - 

scheme can readily be extended to three dimensions via the reference plane approach. , 
This approach involves the definition of a reference plane system in which the threer 
dimensional volume under consideration is spanned by an appropriately selected series • 
of planes which intersect the boundaries of th', considered volume. The equations of 
motion within the reference planes are ex essed in a quasi- streamline coordinate system-,- 
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where quasi- streamlines are the projections of the actual stream surfaces onto these 
reference planes. Then, although the actual streamlines are not traced, streamline 
interpolation procedures are minimized. 

In addition to minimizing steamline interpolation procedures, use of the reference 
plane approach ]^s other distinct advantages. By developing the equations of motion with 
respect to different reference plane systems (Cartesian, cylindrical, and line source), 
complex geometric configurations may be analyzed. In figure 2(a), a reference plane net- 
work is depicted for a typical nozzle module, wherein the line source system shown alle- 
viates the need for adding reference planes as the sidewall opens. The addition or dele- 
tion of reference planes is provided for automatically, based on their proximity to walls. 

A more complex situation is depicted in figure 2(b) for the flow field downstream of the 
modules. For this calculation, a combination of several systems is employed and provi- 
sions are included for automatic switching from one system to another as the character 
of the boundary surfaces changes. The reference plane system also caters to the usage 
of reference plane characteristics at all boundary points. This approach’ is generally 
recognized as the most accurate boundary calculational procedure (ref. 4). However, it 
proves cumbersome when employed in conjunction with nonreference plane networks due 
to the complex interpolation procedures then required. 

The reference plane characteristic technique has been widely used for the calcula- 
tion of three-dimensional supersonic flow fields, and the authors had previously developed 
a program employing this technique for the calculation of nozzle exhaust flow fields (refs. 5 
and 6), which is in current usage at NASA Langley Research Center (refs. 7 and 8). That 
program, as well as most reference plane characteristic (ref char) codes in common usage 
(refs. 9 and 10), employs an inverse scheme wherein interpolations are performed to obtain 
data at the intersection of the quasi- characteristics with the initial data surface. Com- 
parisons of such refchar codes with shock capturing finite difference codes (ref. 11) have 
led to the general conclusion that such difference codes are better able to analyze complex 
flow fields with multiple secondary shocks. From experience gained with the authors' 
original refchar code, it was felt that the Inability to successfully, analyze such flow fields 
was primarily due to the Inverse Interpolation procedures employed. Such procedures 
tend to ignore the presence of weak waves by allowing the quasi- characteristic lines to 
arbitrarily cross each other. The numerical diffusion associated with these interpolations 
can become significant, particularly when the local Courant number (ratio of overall 
marching step to local maximum allowable marching step) is much less than one. The 
smearing of these weak waves is enhanced by resorting to higher order interpolations on 
the initial data line. 

To treat complex multiwave flow fields and still retain the advantages that reference 
plane methods afford, two new numerical codes have been developed. Program CHAR3D 
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is a ref char code which employs the wave preserving network depicted in figurje 3(a) as 
compared to the standard inverse network of figure 3(b). This new network tends to pre- 
serve wave systems and secondary shock waves have been successfully captured with it q 
with a minimum of smearing. In addition, CHAR3D employs a nonisentropic pressure- 
density relation along, streamlines, to- calculate shock entropy losses and. utilizes conser- 
vation variables in constructing derivatives normal to the reference plane. Program 
BIGMAC is a reference plane finite difference code utilizing a quasi- streamline grid in 55 
the reference planes as depicted in figure 3(c). BIGMAC captures shock waves via the 
use of conservation variables in conjunction with a one-sided difference algorithm. 
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SYMBOLS 

equilibrium sound speed, ft/sec (m/sec) 
specific heat at constant volume 

conservation variables (k = 1 to 6 ) defined in text (see eq. ( 1 )) 

conservation variables (k = 1 to 6 ) defined in text (see eq. ( 1 )) 

conservation variables (k = 1 to 6 ) defined in text (see, eq. ( 1 )) 

conservation variables (k = 1 to 6 ) defined in text (see eq. ( 1 )) 

stagnation enthalpy, ft^/sec^ (m^/sec^) 
static enthalpy, ft^/sec^ (m^/sec^) 
defined in text (see eq. ( 1 )) 
index of data point in reference plane 
index of reference plane 
defined in text (see eq. ( 1 )) 
index of marching step 
Mach number in reference plane <lAe 
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ft unit normal to surface 

P ' . pressure, lb/ft2 (N/m2) 

<1 ' ' 'magnitude of velocity in reference plane, ft/sec (m/sec) 

S entropy, ft Vsec^-°R {m^/aec^-K) 

T temperature, ®R (K) 

^ flow velocity vector 

u velocity component in marching direction in reference plane, ft/sec (m/sec) 

V velocity component normal to reference plane, ft/sec (m/sec) 

W velocity component in reference plane normal to marching direction,, 

ft/sec (m/sec) 

Xjr,2 Cartesian coordinates 

tfO,z line source coordinates 

x,0,r cylindrical coordinates 

r equilibrium isen tropic exponent 

p density, slugs/ft^ O^g/m®) 

# fuel-air equivalence ratio 

^ velocity direction in reference jdane, rad 

velocity direction with respect to reference planes, rad 

Arrows over symbols denote vector quantities. Coordinate subscripts denote dif- 
ferentiation with respect to the coordinate. 



GOVERNING EQUATIONS 


Program BIGMAC 

The equations of motion for the steady inviscid flow of a gas mixture in chemical 
equilibrium, written in conservation form with respect to the streamline reference plane 
system described, are: 


Ex + Fy + G^ + H - §2 - F 2 tan a = 0 

where tancK is at constant x and for k= 1 to 6 


( 1 ) 
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Program CHAR3D 

In nonconservation form, these equations (in continuous regions of the flow field) 
may be written 
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( 2 ) 


V • (pV) = 0 ^ . 

p(V • + VP = 0 ^ 

V* VH = 0 

^ = 0 > ^ . 

The eqiiations may be cast in characteristic form with respect to the reference plane 
systems described by writing the above equations in scalar form and transposing those 
terms not involving variations in the x- or z-direction onto the right side. Then, the 
left side is identical to the corresponding two-dimensional system in the X,Z plane. 
The equations in reference plane characteristic form (ref. 6) may be written: 


Along 


_ dz _ cos 0 sin 0 ± ^ 


cos^ i<t>) - 1 
where = (u^ + w^ya| and )3^ = - i 

d0 ± ———K d(ln P) — dx 
TM^ , 


where 


and 




F* = (sin <|> - cos <p) (tan xt^)y + — p 


tan 


(In Pi 


<py tan (cos (^ + X* sin <p) + (j2 “ ^*Jl) tan^ xf/ 


dx = + Jj d(ln x) 

2 ^ 


(3) 


(4) 


d(tan tf^) = jP d(ln P) + G dx 
FM^ 


(5) 

( 6 ) 
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Along 



where 


G = 


-1 

(In P)y 

cos 0 

m2 


+ tan if/ (tan + tan 4^ (l + tan^ cos <p + ‘32 sin 0] 


The flow deflection angles 0 and 0, velocity components u, v, and w, and stream- 
line and characteristic orientations A* are shown in figure 4. A detailed description 
of orientations with respect to the various reference plane systems may. be found in ref- 
erences 5 and 6. ' 

CALCULATIONAL PROCEDURES .! 7 

Interior Point Calculation ' ' 

Properties are desired at the grid point (I,J,K) shown in figure 4. The allowable 
step size Ax is determined by satisfying the CFL condition. For BIGMAC, this requires 
that the intersection of the Mach cone from (f,J,K) with the initial data surface falls within 
the numerical domain as depicted (i.e., the quadrilateral (I,J+1), (I-1,J), (I,J-1), (I+1,J)). 
Note that the effective numerical domain for the characteristic calculation includes the 
points I + 1 and I - 1 on planes J - 1 and J + 1; hence, a larger step may be taken 
with CHAR3D ~ ; ' ' . ^ 

BIGMAC.- The MacCormack (ref. 12) scheme, used to difference equation (1), yields 


%,J = Ei,J - 2 ^(Fy + G^ - F^ tan « + ^ 


(7a) 


where 


tan a = ±2 
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^yi + ^Yo 
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^ ^I,J ' 
~ %J " 


^^2 "" ^I,J+1 " ^I,J 
^^2 = ^U1,J - H,J 


for any variable f and 

E,_J V E, j - 2 to[fy * Gi - ^ F, tan 3^^ 


Where 


(7b) 


tan a * t2 


^I,JT1 - =^I,J \ 
Ayi + Ay^ 




tiit. 



2* T = Zt T + 1 T— ^ ^ 


a,J " ‘^I,J \ho u 


I,J 


*r,j “ “i,j 2 




hgU 


'I,J 


ihgU 


I,J 


Ax 


The physical variables are obtained by the following iterative procedure. A value 
of u is assumed. Then, 

p = E(l)/u (8a) 

P = E(2)-E(l)u (8b) 
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V = E(3)/E(l) 

t - 

(8c) 

• w = E(4)/E(l) • . 

(8d); 

H = E(5)/E(l) 

(8e) 

♦ = E(6)/E(l) 

(8f) 

h = H - •i(u^ + v^ + w^) 

(8g) 

The following three parameter curve fits (based on data from ref. 13) are incorporated ; 
into this code and are described in detail in the appendix of reference 6. 

h = h(P,$,T) 

(9a) 

p = p(P,#,T) 

(9b) 

r= r(p,$,T) 

(9c) 


The value of h obtained in equation (8g) yields T via an inversion of equation (9a). 
Equation (9b) yields an alternate value of the density compared to that obtained in equa- 
tion 8(a). The value of u is perturbed and the procedure repeated until the two values 
of density agree to within a specified tolerance. 

CHAR3D .- Point f in figure 5 is located along the quasi- streamline by the relation 

^f,j = V" V 

where a = 1, b = 0 in the predictor step and ^ ^ in the corrector step. In 

this new wave preserving network, the calculation proceeds upward from the lower bound- 
ary where points (I-1,J) are calculated for ^ reference planes' J to second order prior 
to calculating points (I,J). In addition to the standard initial data array (the points (I,J)), 
an extra array (I,J) is required. To calculate properties at (T,J), the standard initial data 
grid in the reference plane (I-1,J), (I,J), and (I+1,J) is employed to calculate the forcing 
fvinction terms involving derivatives normal to the reference planes. Properties are 
known at points Hj^, Gj, and 1-1 from the calculation of point (I-1,J) to second order. 

Point A is located between and Gj -on the quasi- characteristic X.'^(AI) 
where is defined in equation (3). /dl properties (including forcing functions) are 
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obtained via linear interpolation ^ between Hj and Gj. Then, H 2 is located between 
i and I -I- 1 such that the downrunning quasi-characteristic from G 2 (or B) passes 
through (f,J). To first order, properties at d,J) are calculated using points B and A, 
^ Where Pg and 0g are determined using compatibility relations (eq. (4)) along IB 
and H 2 B. 

Then and <p^ are calculated employing the compatibility relations 


<f>A) 


a 


/_£,\ *bU 


\r«i\ 


\rM^ 


In W 




and 




(<l>i - 0g) - 


.bU 






\TM^ 


in = (aF- + bFf) iSjg 




(10) 


Remaining properties are determined at I via the following streamline relations:.. 


(tan = (tan + 




rM^ 


Hj=Hi 


ai^Hy +b^Hy 
cos <p \cos <f> 


Axg- 


(12) 


<t>l= 0j - 


bSs 4* +b£lL£»„ 

\cos <t> yj^ \co8 <i> yy- 


AXjj- 


(13) 


and in continuous regions of the flow 


(p/pr)j.{p/pr)^-[a[^(p/pr)^ 


+ b 


tan W 


cos 4> 




I Ax. 




(14) 


. ^Linear interpolation along a characteristic line calculated to second order is con- 
sistent to second order. (See appendix.) 
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.The flow velocity is obtained via the relatipn ' v 

r ni/2 

- (15) 

where T- is obtained via an Inversion of equation (9b) ^with p_, P-, and known^ 
and h^ is obtained employing equation (9a). Then, Fj^ is obtained from equation (9c) 
and a| = T-Pj^p-. 

This calculation is performed for points I in ^ reference planes to first order. 
Then, cross derivatives a/ ty are evaluated at f employing the relation 


where 








(ay, + ayj) 




(16) 


'">'1 " ■ >'w-l '">'2 = V*l ■ ^i,J 

•an o = te) 

' /x,T7 

Im\ ^f,j " 

' ' " ^I,J “ "“f-M : 

Derivative's are made the same way at the initial station I, except here ' 8f/az is eval- 
uated by ■ ■ • ' 


^4,y "‘i+iyj" ; 

CHAR3D,,in addition to the centered difference algorithm described above, has the 
option of evaluating, cross derivatives via an alternating one-sided difference algorithm. 
For this option, derivatives are evaluated as described in the section for BIGMAC. Cross 
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derivatives are required for the variables ■ P, <p, H, and In evaluating 

cross derivatives for P, <l>, and \{/, conservation variables are employed as follows: 

•Py = F(3)y- F(l)Vy- VF(l)y . 

' Wv cos <|) - Uy sin 0 I • • . . ■ . 

qvy - (uuy + wwy) tan yp 

■ (tani^)y = 

^ q2. - J 

where ^ - 

u = q cos <p w = q sin <f> v = q tan \f/ 

q2 - u2 + 

and ; • 

E(3)y-vE(l) ^ 

V" W) , 

E(2)y - hghgPy - U E(l)y 

i^) - ; 

E(4)y - w E(l)y. 

■ E(i) 

The conservation variables E(k) and F(k) are given by equation (1). The use of con- 
servation variables in construction of these cross derivatives tends to suppress oscilla- 
tions, that occur when employing physical variables to difference across shock waves. . 
However, the use of a one-sided difference algorithm in conjunction with CHAR3D. tends 
to produce spurious results in regions of large cross flow. 

In the characteristic reference plane algorithm, cross flow variations are expressed 
via the forcing function terms F* appearing in the right side of the compatibility rela- 
tions (eq. (4)). These terms are assumed to vary mildly within an integration step. When 
a one-sided algorithm is employed to evaluate cross derivatives in the vicinity of shocks, 
the values of the forcing function terms may vary greatly between the predictor and cor- 
rector stei>s. In addition, the numerical domain of dependence is somewhat vague for the 
characteristic reference plane approach in conjunction with one-sided differences, so that 
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part of the problem may be due to s^bility. . The recommended approach for evaluating 
cross derivatives in CHAR30 is to employ conservation variables in conjunction with a 
centered difference algorithm, although thi^ matter requires further study. 


In CHAR3D, secondary shocks are captured as rapid changes spread over approxi- . 
mately three grid points. These discontinuities are preserved by use of the wave network 
described which performs all interpolations off quasi- characteristic surfaces. The 
entropy change associated with these shocks is evaluated employing a nonisentropic 
pressure-density relation (illustrated here for a perfect gas) 

V- Vin(p/pr)= V-^ (19) 


For a shock of strength | (pressure ratio; across shock), this change is deter inined 
employing the relation (for perfect gas) 


Cv 


In I - rin! 


(r+ 1 )^ -«• (r- 1 ) 
_(r- 1)1 + (r + i)_ 


(20) 


where is the entropy change along a streamline produced by the captured shock. 
This relation involves only die pressure distribution in the vicinity of the shock and is 
readily applied in regions of noninteracting shocks as follows. Let 


F(^,D = 


(r -t- i)g -H (r- 1 ) 
(r- 1)4 + (Ft 1) 


Assume a Shock is spread over the marching interval K - 1 to 6 (fig. 6) for a typical 
quasi- streamline. Then 1 represents free stream conditions for this shock. The 
entropy change in the interval K - 1 to K is then expressed by 


m = 7 ^ 


yy 


^VJ 


AS) 


= In 


K-1,K V '1,K ' '14C-1 


-^1,K /^1,K-1 ) 
^1,K-1\ ^1,K 


where 


*1,K ' % 




Then 
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cos <P \ ' ■ /y 


exp 




K-1,K 



Since the shock geometry does not appear in the entropy jump relation, the entropy rise 
associated with extremely complex three-dimensional shocks can be accurately obtained. 
Special provisions have been incorporated into the program for the computation of singu- 
lar points at the juncture of intersecting shock waves and/or shock reflection points. At 
such points the streamline undergoes a discontinuous pressure. rise corresponding to that 
through both shock waves, the shock intensities are different, an entropy discontinuity 
occurs separatii^ the different zones, and a vortex of infinite intensity results. Numeri- 
cally, the entropy procedure, described would predict an entropy rise associated with this 
pressure jump. Theoretically, this occurs in the limit of vanishing mass flow, while 
numerically the finite mass within this region would lead to unduly large entropy levels. 
Special coding has been incorporated at such singular points to suppress these "numerical" 
peaks. 

Wall Point Calculation 

Solid surfaces are prescribed via discrete contour data and fitted via a newly devel- 
oped method based on the use of partial cubic splines (ref. 14). The surface fitting is done 
by a separate geometry package and the array of coefficients generated is stored on tape.. 
BIGMAC and CHAR3D employ this coefficient data in conjunction with a surface interpola- 
tion procedure yielding highly accurate values of the dependent variable and surface unit 
normal. 

In both BIGMAC and CHAR3D, wall point calculations are performed employing a 
reference plane characteristic calculation. In figure 7, CD is the intersection of the 
reference plane y = Vq with the surface z = f(x,y). Reference planes are oriented so 
that the surface normal lies nearly within the reference plane. For sidewall calculations, 
this is accomplished via local coordinate rotations. 

In CHAR3D, and are evaluated utilizing the characteristic compat- 

ibility relation (eq. (4)) along BC, the normal momentum equation (eq. (6)) along the 
streamline projection CD*, and the relation V • n = 0 applied at C, which yield the 
relation 

sin 0c = (fx)c cos 0c + (fy)c tan 0c (21) 

The compatibility equation yields a relation between Pc and <pQ, and the normal 
momentum equation yields a relation between Pc and 0c . This system is solved in 
the context of the wave preserving network previously described by a simple iterative 
procedure. 
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Li BIGMAC, this iterative procedure is eliminated by combining the normal momen* 
turn equation with the quasi* streamline momentum equation, yielding the follo^^g system 
of equations for P^, (w/u)^,and (v/u)^: 


(32) 




(w/u)^. 



‘ ® -Wc 


In Pc 


^2 



(v/u)c 


Rs 

* 

Q 

1 






where 


I 2 \ t \ 

r" - * in Pg * WbC Wb 

^2 * (^x)c 


^3 


u j®SL • tan ^ (B cos 4 > + C sin 0)1 Ax 
pu h^h2h3 L j 




and 




A = hih 


V I ^ I T E(l) , G(l) 
jh, -T + PVy + Jl -j^ ♦ J2 -J- 


B - P(2)y - U F(l)y - J, 


C . F(4)y - W FO)y - Jj 


V F(l) 
1»3 


Egi, . F(3)y - y F(l), W, M ^ Jj ^ ' 

Then, relations applied along the streamline projection CD* yield remaining flow var- 
iables at C, in conjunction with the equilibrium curve fits (eq. (9)), for both programs. 
The process is then repeated with coefficients averaged for second order accuracy. 
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;-;i: ; tv: : Interior 'Corner ■ 

s:rr.'- r ;v-,v. .• . . 

Interior corners occur in the internal modules and are discretely treated as the 
intersection of specified surfaces, as depicted in figure 8. A detailed description of these 
corner calculatipns with respect to the various reference plane coordinate systems may 
be found jih reference 6. The procedure is outlined here for a Cartesian system where 
'fee intersecting surfaces are prescribed by z = f(x,y) and y = g(x,z). 

The relation V • n = 0 applied to both intersecting surfaces at C (fee point to be 
calculated) yields the flow deflection angles . and explicitly^ 



tan"^ 


} ’ SzV/ 






1 " fySz / 


(23) 




Then, a redundant procedure is employed wherein reference plane calculations for fee 
pressure at C are performed in the reference planes z = Zq a,nd y = y^- This yields 
two values of pressure and which differ due to evaluating .the cross deriva- 

tive forcing function terms in the compatibility relations via backward differences. A 
weighting of these pressures is performed by accounting for the relative wave strengfes 
in each of these reference planes. This gives fee stronger weighting to the calculation 
performed in the reference plane containing the dominant waves via the relation 


^Ci 


^ + ^^A2C ^^2 


(24) 


Streamline relations are performed along the corner CD, and the process is repeated for 
second order accuracy. 


Shock Point Calculation 

A discrete three-dimensional shock point calculation is performed for the nozzle 
imderexpansion shock, which propagates into fee nonuniform external stream surrounding 
fee vehicle. In figure 9, subscript 2 ’ refers to the shock free stream. Shock geometry 
is defined in terms of fee direction cosines of a and 8, where ^ is fee angle made by 
fee shock cut wife fee reference plane and d is fee crosscut ah^e. "For given-values 
of O’- and /3, the shock normal is 
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(25) 


fig a lx COS a sin ^ - iy sin a + iz cos a cos ^ 


where ix>iy>iz the^unit vectors in the x-, y- , and z-directions. The characteristic - 
relations on the free stream side yield flow properties at C 2 . The Hugoniot relations 
in a shock normal system yield properties at C^. The compatibility relation along AjCi 
yields an alternate value of pressure Pcj* The angle /3 is perturbed locally until Pq^ 

from the jump relations equals from the compatibility relations to within a speci- 


fied tolerance. This procedure is performed in all reference planes, and the process is 
then repeated using updated values of the crosscut angle a. The complete details of this 
procedure including rotation into the shock oriented system, jump relations, and iterative 
procedures may be found in references 5 and 6. 


Contact Surface Calculation 


A three-dimensional contact surface is significantly more complex than its two- 
dimensional counterpart, since the streamlines on each side of the discontinuity not only 
differ in velocity magnitude but also may be highly skewed with respect to each other. In 
figure 10, <x and /3 are as previously defined for the shock calculation, and the stream- 
lines passing through C emanate from* Dj on the lower side and D 2 on the upper 
side. Hence, discontinuities exist in the flow angles 0 and . ^ at point C. The bound- 
ary relation V ♦ n = 0 applied at Cl and C 2 yields the relations 


sin - ^Ci) ^ tan Q! tan 
sin + tan Q! tan \{/q^ 



(26), 


Then, characteristic compatibility relations may be applied along AiCi and B 2 C 2 
yielding ^Ci “ ^Ci ^^2 ” ^€2 The normal momentum relations 

applied along the streamline projections CiDj^ and C 2 D 2 yield relations between 

and " ^€ 2 ' ^ given value of the crosscut angle a, a value of ^ 

is obtained via an iterative process satisfying the above relations and the boundary condi- 
tion Pqj^ = Pc 2 * This procedure is i)erformed in all reference planes and repeated with 

updated values of the crosscut angle a. Again, complete details may be found in refer- 
ences 5 and 6. • ’ 


RESULTS 

Internal corners represent just one segment of the overall boundary calculational 
procedure and hence must be calculated as part of the overall marchi^ procedure. 
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Recently, inviscid corner flow fields have been studied in detail (refs. 15 and 16) utiliz- 
ing conical coordinates in a timelike marching procedure until conical invariance was 
achieved. While these schemes do yield the flow field devils in the corner region, they 
.are not, applicable to general three-dimensional flow problems which are noncOnical. . 

Corner results are presented using the general interior corner point calculation 
outlined above and previously described in references 5 and 6. Results for a 5° double 
expansion corner are depicted in figures 11 and 12. These results were obtained with 
CHAR3D starting from uniform initial flow conditions = 845, = 2.94, <p s \}/ s Oj 

with an 11 x 11 Cartesian grid. Results are shown after nine axial marching steps and 
the axial pressure variation at the corner is also indicated. Similar results have been 
obtained with BIGMAC. 

An expansion- compression has been calculated using BIGMAC which yielded the 
results depicted in figure 13. These results were obtained with an 11 x 11 Cartesian 
grid for initially uniform flow (Me© = 2) and are depicted after 10 axial marching steps. ■ 
Results are compared with the detailed solution of Shankar (ref. 16) and the experimental 
results of Nangia (ref. 17). 

Results for the double compression corner, as obtained by BIGMAC, are shown in 
figure 14 after 35 axial marching steps. A 12 x 12 line source network was employed 
with initially uniform flow at Mo© = 3.17. A comparison is made with Shankar's numer- 
ical results (ref. 16) and the experimental results of Charwat and Redekeopp (ref. 18). 

The above results verify the accuracy and validity of the interior corner procedure 
employed and, hence, yield credibility to the application of this procedure for general cor- 
ner calculations within "truly" three-dimensional flow fields. 

To demonstrate results obtainable with the new wave preserving network of CHAR3D, 
a simple two-dimensional inlet flow field is calculated. Calculation was performed with 
a uniform equally spaced initial profile (P « 845, M » 2.94) employing 11 and 21 grid 
points. Wall pressures are depicted in figures 15 and 16 for three shock reflections. 

After the fourth reflection, the flow on the upper boundary is subsonic, and thus the pro- 
gram could not calculate past this region. Note that both the pressures obtained as well 
as the propagation rates are in excellent agreement with the exact solution and no addi- 
tional smearing results from wall reflections. 

A complex internal module flow field calculation (square nozzle) as depicted in fig- 
ure 17 has been performed using BIGMAC. This flow field is characterized by the initial 
interactions of expansion waves emanating from mutually perpendicular surfaces and the 
subsequent interaction of enveloping shock systems generated by recompression on the 
upper wall and sidewall. This calculation employed a 21 x 11 Cartesian network, with 
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additio^l. reference planes being inserted as, the. sidewall opened. At the straight section 
the final network was 21 x 18. Uniform flow properties = 845, .= 2. 94,1 were 

stipulated at the nozzle entrance. Pressure contours pn the symmetry plane are dei|icted 
in figure 18. Of particular interest is the intersection of four three- dimensional, shpck ^ 
surfaces.at x » 17 and z ,= y = 0. This results from the reflection of the envelope, 
shock produced bye the sidewall and the r.eflection of the envelope shock produced by .the 
upper wall, resulting in an approximate 15/1 pressure ratio at this location. The axial 
pressure variation along the corner is depicted in figure 19 and pressure variations along 
several streamlines in the symmetry plane are depicted in figure 20. 

All results presented employed a perfect gas option with .F = 1.4 for the sake of 
simplicity. . 'The equilibrium option has been extensively used and tested (refs. 5, 6, 7> 
and 8) and provides.np further insight into these problems. The results wpre all obtained 
with relatively crude grid networks, yet provided accurate and detailed flow field results. 
Further grid refinement would yield somewhat better flow resolution, if desired or neces- 
sary. It should be noted that due to the use of disc storage techniques, as employed in 
both programs, flow field resolution is riot limited by machine core storage. > . • 

CONCLUDING REMARKS . ‘ • 

Two new computer codes have been developed for analyzing complex three- * 
dimensional supersonic flow fields. Their use of a quasi- streamline network in conjunc- 
tion with a reference plane grid allows for the calculation of complex geometric config- 
urations and caters to highly rotational, variable composition flow fields. Both BlGMAC 
and CHAR3D are currently running internal flow codes with perfect gas or equilibrium 
hydrogen-air chemistry options. 

CHAR3D employs a totally new grid network which caters to both the following of 
streamlines and the preservation of wave systems. This is done in conjunction with an 
axial marching procedure. Hence, in addition to its application to three-dimensional ref- 
erence plane systems, it is equally applicable to "viscous" characteristic techniques, 
since forcing functions are also employed. 

BlGMAC employs the commonly used MacCormack algorithm in conjunction with 
conservation variables and hence falls in the general classification of finite difference 
shock capturing codes. However, it does this in conjunction with a reference plane 
streamline grid which provides significant advantages for the flow fields treated. 

Both programs treat complex three-dimensional flow fields accurately, locating 
secondary shock waves and evaluating flow field properties in their vicinity Including 
wall and interior flow entropy. From our limited experience with these codes, CHAR3D 
appears best suited to flow fields wherein the predominant wave propagation occurs within 
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the reference planes. For such flow fields,' CHAR3D Mth half the ^id points yields - 
results' comparable to those of BIGMAC. In addition, no overshoots occur in the vicinity' 
of "shock waves and a larger marching step may be taken. However, when the assumption 
of a Mildly varying forcing function is violated (i.e., in the vicinity of strong crosswise 
compressions) BIGMAC would be the preferred code. This program has no preferen- 
tial direction and has been shown capable of calculating arbitrary multishocked three- 
dimensional flow fields. 

Our current effort is devoted to extending both these codes for the calculation of 
the flow field downstream of the engine modules. This calculation is performed in the 
authors' previous code and similar procedures will be incorporated.' 'Future efforts will . 
involve the incorporation of finite-rate hydrogen-air chemistry, frozen chemistry, and . 
associated sudden freezing criterion. In addition^ the extension of these codes to mixing 
calculations along the plume interface is anticipated. '■ 

It should be noted that while the calculation of nozzle exhaust flow fields has been 
specifically discussed, both codes are capable of analyzing quite general three-dimensional 
flow fields. Results to date, indicate that these techniques yield minimum smearing of cap- 
tured shocks, even after multiple reflections and/or intersections. Thus, these cpdes 
appear capable of calculating inlet type flow fields and can readily be modified, to calculate 
the simpler problem of external supersonic flows. 



APPENDK 

LINEAR INTERPOLATION ON CHARACTERISTIC " ' “ 

A SECOND ORDER PROCEDURE 

In previous reference plane characteristic codes (employing inverse interpolation 
procedures), data are interpolated on a noncharacteristic surface. To achieve full second 
order accuracy, most codes resort to higher order interpolation procedures. Such pro- 
cedures are helpful in smooth regions of the flow field but are detrimental in regions of 
weak discontinuities. In such regions, linear interpolation is more accurate as explicitly 
discussed by Sedney (ref. 1). The authors had performed an independent study (unpub- 
lished) on such higher order interpolation procedures and concluded that for general 
miiltiwave flow fields, linear interpolations provide the most accurate results. 

Now, with this new ’*wave preserving” network, all interpolations are performed on 
characteristic lines. Employing a linear interpolation procedure on a characteristic line 
calculated to second order is consistent with a second order algorithm. This point can 
be inferred from Ferri's article (ref. 19) but apparently , is not universally accepted. (See 
reL 1.) Hence, a simple proof of this statement is presented. 

Along any line AC, a series expansion for the pressure and flow deflection are 
written tr ■ .« 


Pc = Pa + (Px)a + (Pxx)a 2 + 0 (Ax)3 

(Al) 

^C = <^A + Wa ^ + (^xx)a (^')V2 + 0 (Ax)3 

(A2) 

but. ^ - 


(^x)c “ (Px)a (^xx)a ® ■ 

' (A3) 

(''■x)c = (^x)a ^ (■'■xxV ^ 

(A4) 

where x denotes distance along AC. 

.. A 
4 ' *■ 

Substituting -equation (A3) into equation (Al) and equation (A4) into equation (A2)~ 

results in 


Pc = Pa + \{^ x)a + (Px)cl + 9 

(A5) 

0c = <#>A + \{^ x)a + + 0 (Ax)3 

(A6) 
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APPENDIX — Continued 


The previous expressions are valid along ahy line AC. Assuming that AC is a down- 
running characteristic, the compatibility relations at points A and C are 


Wa - *a(Px)a = « 
■ K)c -A:(Px)c = 0 


Where.; 


■'•A = 


/3 


FM^P 


Solving the system of equations (A5) to (A8).for (Px)a (^x)c resiilts in 

. (Px)a.° Aa) |»C - »a) - Ac (Pc - Pa)] 


and 


(Pxfc = ax(Ac-Aa) ['fc - * a ] - Aa(Pc - Pa) 


(A7) 

(A8) 


(A9) 


(AlO) 


Now consider a point x* between A and C. The pressure at this point to sec- 
ond order is given by 

P* = PA*(Px)A(*’-*A)-^^^T3^;^(**-*Af (All) 


where (Px)a (^x)c given by equations (A9) and (AlO). Up to this point, all 

relations are quite general and have not required that a second order compatibility rela- 
tion exist between A and C. We now make use of this relation by stating that by sec- 
ond order, we imply that the relation 

(♦c - •Ca) - (^^T^)(Pc - Pa) = » <A12) 

is satisfied between points A and C in a convergent fashion as detailed in reference 19. 
Then, substituting equation (A12) into equations (A9) and (AlO) results in 


(Px)c - (Px)a = » 


(A13) 
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and 


APPENDIX - Concluded 

V 



Pc -P a 

*C-*A 


Hence, substituting these relations into equation (All) yields 


V- (A14) 


which clearly demonstrates that a linear interpolation for pressure (or flow deflection) 
on a characteristic Adulated to second order is consistent with a fully second order . 
approach. 


•P’ - Pa * ( 


* « (Ax)® 
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(a) lii line’ 'source system for internal nozzle module. JW designates sidewall. 











□ WAVE GR® POINT 
O STREAMLINE GRID POINT 

^ ~ ' r 


(a) Reference plane grid network for CHAR3D. 




(b) Standard reference plane (c) Finite difference network, 

characteristic network. 

Figure 3;- Reference plane networks. 
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REFERENCE PLANE J (Y >CONSTANT) 



O STREAMLINE GRID ARRAY 

• INTERMEDIATE POINT 

(CHARACTERISTIC 

CALCULATION) 

□ EXTRA CHARACTERISTIC 
ARRAY (INTERPOLATED 
ALONG CHARACTERISTIC) 




Figure 5.- CHAR3D interior point grid. 
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Figure 9.* Shock surface calculation. 


Figure 10.- Contact surface calculation. 
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Figure 12.~, Pressure distribution for 5® expansion corner. PM^ xneans Prandtl- Meyer. 
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Figiire 14.« Compression corner. - 6^ is wec^e angle; Pj is total pressure. 
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Figure 19.- Streamline pressure distribution at sidewall corner of square nozzle. 

21x11 GRID M<jo««2.94 




